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Chapter 2 

Theoretical Models of SOC 
Systems 

by Markus J. Aschwanden 

How can the universe start with a few types of elementary particles at the big bang, 
and end up with life, history, economics, and literature? The question is screaming out 
to be answered but it is seldom even asked. Why did the big bang not form a simple 
gas of particles, or condense into one big crystal? (Bak 1996). The answer to this fun- 
damental question lies in the tendency of the universal evolution towards complexity, 
which is a property of many nonlinear energy dissipation processes. Dissipative non- 
linear systems generally have a source of free energy, which can be partially dissipated 
whenever an instability occurs. This triggers an avalanche -like energy dissipation event 
above some threshold level. Such nonlinear processes are observed in astrophysics, 
magnetospheric physics, geophysics, material sciences, physical laboratories, human 
activities (stock market, city sizes, internet, brain activity), and in natural hazards and 
catastrophes (earthquakes, snow avalanches, forest fires). A tentative list of SOC phe- 
nomena with the relevant sources of free energy, the physical driver mechanisms, and 
instabilities that trigger a SOC event are listed in Table 2. 1. 

A prominent theory that explains such nonlinear energy dissipation events is the 
so-called Self- organized criticality (SOC) concept, first pioneered by Bak et al. (1987, 
1988) and simulated with cellular automaton models, which mimic next-neighbour in- 
teractions leading to complex patterns. The topic of SOC is reviewed in recent reviews, 
textbooks, and monographs (e.g., Bak 1996; Jensen 1998; Turcotte 1999; Charbon- 
neau et al. 2001; Hergarten 2002; Sornette 2004; Aschwanden 2011a; Crosby 2011; 
Pruessner 2012). SOC can be considered as a basic physics phenomenon - univer- 
sally occurring in systems with many coupled degrees of freedom in the limit of in- 
finitesimal external forcing. This theory assumes a critical state that is robust in the 
sense that it is self-organizing, like a critical slope of a sandpile is maintained under 
the steady (but random) dropping of new sand grains on top of the pile. Individual 
avalanches occur with unpredictable sizes, uncorrelated to the disturbances produced 
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Table 2.1: Examples of physical processes with SOC behavior. 



SOC Phenomenon 


Source of free energy 


Instability or 




or physical mechanism 


trigger of SOC event 


Galaxy formation 


gravity, rotation 


density fluctuations 


Star formation 


gravity, rotation 


gravitational collapse 


Blazars 


gravity, magnetic field 


relativistic jets 


Soft gamma ray repeaters 


magnetic field 


star crust fractures 


Pulsar glitches 


rotation 


Magnus force 


Blackhole objects 


gravity, rotation 


accretion disk instability 


Cosmic rays 


magnetic field, shocks 


particle acceleration 


Solar/stellar dynamo 


magnetofriction in tachocUne 


magnetic buoyancy 


Solar/stellar flares 


magnetic stressing 


magnetic reconnection 


Nuclear burning 


atomic energy 


chain reaction 


Saturn rings 


kinetic energy 


collisions 


Asteroid belt 


kinetic energy 


colhsions 


Lunar craters 


lunar gravity 


meteroid impact 


Magnetospheric substorms 


electric currents, solar wind 


magnetic reconnection 


Earthquakes 


continental drift 


tectonic shpping 


Snow avalanches 


gravity 


temperature increase 


Sandpile avalanches 


gravity 


super-critical slope 


Forest fire 


heat capacity of wood 


lightening, campfire 


Lightening 


electrostatic potential 


discharge 


Traffic collisions 


kinetic energy of cars 


driver distraction, ice 


Stockmarket crash 


economic capital, profit 


political event, speculation 


Lottery win 


optimistic buyers 


random drawing system 



by the input. Sandpile avalanches are a paradigm of the SOC theory, which has the fol- 
lowing characteristics: (1) Individual events are statistically independent, spatially and 
temporally (leading to random waiting time distributions); (2) The size or occurrence 
frequency distribution is scale-free and can be characterized by a powerlaw function 
over some size range; (3) The detailed spatial and temporal evolution is complex and 
involves a fractal geometry and stochastically fluctuating (intermittent) time character- 
istics (sometimes modeled with 1/f-noise, white, pink, red, or black noise). 

There are some related physical processes that share some of these characteristics, 
and thus are difficult to discriminate from a SOC process, such as turbulence, Brownian 
motion, percolation, or chaotic systems (Fig. 2.1, right column). A universal SOC 
theory that makes quantitative predictions of the powerlaw-like occurrence frequency 
and waiting time distributions is still lacking. We expect that the analysis of large 
new databases of space observations and geophysics records (available over at least a 
half century now) will provide unprecedented statistics of SOC observables (Fig. 2.1, 
left column), which will constrain the observed and theoretically predicted statistical 
distribution functions (Fig. 2.2, middle column), and this way confirm or disprove the 
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Observables 



Statistical Distributions 



Pliysical Processes 



Photon counts 
Fluence 
Duration 
Time intervai 
Lengtii scaie / 
Area \ 
Volume 
Energy 

Current density 
Magnetic fieid 
Emission measure 
Density 
Temperature 



:> 



Frequency distribution 
Waiting time distribution 



White noise spectrum 
Pinl< noise spectrum 
Brown noise spectrum 
Biacl< noise spectrum 



Poweriaw distribution 
Exponential distribution 
Log-normal distribution 
Binomial distribution 
Pareto distribution 



Poisson distribution 



Self-organization 
Forced Criticality 
Brownian Motion 
Tsailis Entropy 
Hyper-diffusion 
Levy Flights 



Stationary SOC 
Nonstationary SOC 
Fractal-diffusive SOC 



Phase Transisions 
MHD turbulence 
Percolation 



Chaotic systems 



Figure 2. 1 : Metrics of observables, statistical distributions, and physical models that 
need to be defined in order to discriminate SOC from non-SOC processes. 

theoretical predictions of existing SOC theories and models (Fig. 2.2, right column). 
Ultimately we expect to find a set of observables that allows us to discriminate among 
different SOC models as well as against other nonlinear dissipative processes (e.g., 
MHD turbulence). In the following we will describe existing SOC models and SOC- 
related processes from a theoretical point of view. 



The theoretical models can be subdivided into: (i) a mathematical part that deals with 
the statistical aspects of complexity, which is universal to all SOC phenomena and 
essentially "physics-free" (such as the powerlaw-like distribution functions), and (ii) a 
physical part that links the observable parameters to a particular physical mechanism 
(for instance in form of scaling law relationships between physical parameters). 

2.1.1 Statistical Aspects 

The original concept of self-organized criticality was pioneered by Bak et al. (1987, 
1988), who used the paradigm of a sandpile with a critical slope to qualitatively illus- 
trate the principle of self-organized criticality. The first theoretical description of the 
SOC behavior of a nonlinear system was then studied with a large number of coupled 
pendulums and with numerical experiments that have been dubbed cellular automaton 
models (CA-SOC). In essence, a SOC system has a large number of elements with 
coupled degrees of freedom, where a random disturbance causes a complex dynamical 
spatio-temporal pattern. Although such a nonlinear system obeys classical or quantum- 
mechanical physics, it cannot be described analytically because of the large number of 
degrees of freedom. The situation is similar to the N-body problem, which cannot be 
solved analytically for more complex configurations than a two-body system. However, 
such SOC systems can easily be simulated with a numerical computer program, which 
starts with an initial condition defined in a lattice grid, and then iteratively updates the 
dynamical state of each system node, and this way mimics the dynamical evolution 
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and statistical distribution functions of SOC parameters. Hence, such numerical lattice 
simulations of SOC behavior (i.e., cellular automatons) were the first viable method to 
study the SOC phenomenon on a theoretical basis. 

What are the ingredients and physical parameters of a cellular automaton model? 
What is common to all cellular automaton models of the type of Bak, Tang, and 
Wiesenfeld (1987, 1988), briefly called BTW models, is: (1) A S-dimensional rectan- 
gular lattice grid, say a coordinate system .Ti J jt, i = l,...,n:j = l,...,n, fc = 
for a 3-dimensional case S = 3 with grid size n; (2) a place-holder for a physical quan- 
tity Zij^k associated with each cellular node Xij^k, (3) a definition of a critical threshold 
Zcrit, (4) a random input Azij,k in space and time; (5) a mathematical re-distribution 
rule that is applied when a local physical quantity exceeds the critical threshold value, 
for instance a critical slope of a sandpile, which adjusts the state of the next-neighbor 
ceUs, i.e., 

Zi,j,k = Zi,j,k + 1 initial input 

Zi,j,k = Zi,j,k - ^ if^i,j,fe>8, (2.1) 

Zi±l,j±l,k±l = ^i±l,j±l,fe±l + 1 

for the 8 next neighbors in a 3D-lattice grid, and (6) iterative steps in time to update 
the system state Zij^^ {t) as a function of time t. Such cellular automaton simulations 
usually start with a stable initial condition (e.g., an empty system with Zi_j^k{t = 0) = 
0), and have then to run several million time steps before they reach a critical state. 
Once they reach the critical state, avalanches of arbitrary sizes can be triggered by the 
random input of an infinitesimal disturbance Az^ j ^(t), separated by quiescent time 
intervals in between, in the case of slow driving. Statistics of the occurrence frequency 
distributions of avalanche sizes or durations yields then the ubiquitous powerlaw-like 
distribution functions. Fig. 2.2 shows an example of a 2D-lattice simulation, displaying 
a snapshot of a large avalanche (top left panel) as well as the statistical distribution 
function of avalanche sizes L (left) and durations T (right panels). 

Thus, such SOC cellular automaton models primarily quantify the spatio-temporal 
dynamics of complex systems and the related statistics of emerging spatio-temporal 
patterns. They can predict the distribution function N{L)dL of the size or length scale 
L of dynamical events (also called "avalanches"), or the time scale T of an avalanche 
event. Secondary dynamical parameters can be derived, such as the instantaneous 
avalanche volume V{t) or its change dV{t) / dt during an avalanche event, or the total 
time-integrated volume of the avalanche V(t)dt. However, such SOC cellular au- 
tomaton models are "physics-free", because the dynamics of an event is described by a 
mathematical rule that substitutes for the unknown physical process, but is thought to 
be an universal property of SOC systems in a critical state. Applying cellular automaton 
simulations to observed phenomena requires than to substitute a physical quantity into 
the "place-holder" variable Zij^k, which is otherwise undefined in cellular automata, 
except for defining the instability criterion Zij^k > Zcrit in terms of a critical threshold 

Zcrit ■ 

There is a large industry of cellular automaton models, mostly modified versions of 
the original BTW-model, such as the lattice-gas model (Jensen 1998), Conway's game 
of life model (Bak, Chen and Creutz 1989), traffic jam simulations (Nagel and Paczuski 
1995), multiple- strategy agent-based models applied to the financial market (Feigen- 
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Figure 2.2: Examples of a fragmented avalanche (top left) occurring in a 2-D (computer) sand- 
pile (top right) and occurrence frequency distribution of avalanche cluster sizes (left panels) 
and avalanche durations (right panels) of the original BTW sandpile cellular automaton simula- 
tion. The simulations have been performed for a 50 x 50 2-D lattice (middle panels) and for a 
20x20x20 3-D lattice grid (bottom panels). The powerlaw slopes are as = 1.0 and ot = 0.42 
for the 2-D grid (middle panels) and as = 1-37 and ar = 0.92 for the 3-D grid. Reprinted 
from Bak, Tang, and Wiesenfeld (1987, 1988) with permission; Copyright by American Physical 
Society. 
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baum 2003), punctuated equilibrium models applied to biophysics (Bak and Snep- 
pen 1993), sUder-block spring models applied in geophysics (Burridge and Knopoff 
1967), forest-fire models (Malamud et al. 1998; reviewed in Turcotte 1999), applica- 
tions to magnetospheric substorms (Takalo et al. 1993; 1999a,b), to solar flares (Lu and 
Hamilton 1991; Charbonneau et al. 2001), and to stellar accretion disks (Mineshige et 
al. 1994a,b; Pavlidou et al. 2001). More complete reviews of such cellular automaton 
models are given in Turcotte (1999) and Aschwanden (201 la. Section 2). 

2.1.2 Physical Aspects 

Every theoretical model needs verification by experiments and observations. Since a 
cellular automaton model is a purely mathematical model of complexity, similar to the 
mathematical definition of a fractal dimension, its predictions of a particular mathe- 
matical distribution function (such as a powerlaw) may represents a universal property 
of SOC processes, but does not represent a complete physical model per se that can be 
used for applications to real observations. The physical aspect of a SOC model is hid- 
den in the "place-holder" variable Zij^k and its scaling laws with observables, which 
require a specific physical mechanism for each observed phenomenon (for examples 
see second column in Table 2.1). 

If we go back to the original SOC concept of avalanches in a sandpile, say in the ID- 
version, the place-holder variable Zi has been interpreted as a vertical altitude or height 
difference Zi = h i+i — hi between next-neighbor nodes, so that the instability threshold 
corresponds to a critical slope Zi > Zcrit = {dh/dx)crit- For the 3D version, the re- 
distribution rule is given in Eq. (2.1). The instability threshold essentially corresponds 
to a critical point where the gravitational force exceeds the frictional force of a sand 
grain. 

In appUcations to solar flares or magnetospheric substorms , the "place-holder" vari- 
able Zi,j_k has been related to the magnetic field Bijj- = B(xi,j_k) at location Xij,k, 
and a related magnetic energy Eb = B"^ / (Stt) can be defined. The instability criterion 
AS > Zcrit involves than a magnetic gradient or magnetic field curvature (Charbon- 
neau et al. 2001), 



where the sunmiation includes all 2 x 5 nearest neighbors ( "nn") in a Cartesian S- 
dimensional lattice. 

More realistically, physics-based models have been attempted by applying the magneto- 
hydrodynamic (MHD) equations to the lattice field i3(x), obeying Ampere's law for 
the current density j. 




nn=l 




(2.2) 




(2.3) 



which yields together with Ohm's law the induction equation. 



— - = V X (v X B) -I- rjV^B , 



(2.4) 
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and fulfills the divergence-free condition for the magnetic field, 

V • B = . (2.5) 

The instability threshold can then be expressed in terms of a critical resisitivity rj. Such 
cellular models with discretized magneto-hydrodynamics (MHD) have been applied 
to magnetospheric phenomena, triggering magnetospheric substorms by perturbations 
in the solar wind (Takalo et al. 1993; 1999), as well as to solar flares (VassiUadis et 
al. 1998; Isliker et al. 1998). Initially, solar flares were modeled with isotropic mag- 
netic field cellular automaton models (Lu and Hamilton 1991). However, since the 
plasma-^ parameter is generally smaller than unity in the solar corona, particle and 
plasma flows are guided by the magnetic field, and thus SOC avalanches are more suit- 
ably represented with anisotropic ID transport (Vlahos et al. 1995). Other incarnations 
of cellular SOC models mimic the magnetic field braiding that is thought to contribute 
to coronal heating (Morales and Charbonneau 2008). On the other side, SOC mod- 
els that involve discretized ideal MHD equations have been criticized to be inadequate 
to describe the highly resistive and turbulent evolution of magnetic reconnection pro- 
cesses, which are believed to be the driver of solar flares. 

Besides the magnetic field approach of the a physical variable in SOC lattice sim- 
ulations, which predicts the magnetic energy Eb oc B'^/{8tt) or current j oc V x B 
in each node point, we still have to model the physical units of observables. The mag- 
netic field or electric currents can only be observed by in-situ measurements, which 
is possible for magnetospheric or heliospheric phenomena with spacecraft, but often 
observables can only be obtained by remote-sensing measurements in form of photon 
fluxes emitted in various wavelength ranges, such as for solar or astrophysical sources. 
This requires physical modeling of the radiation process, a component that has been 
neglected in most previous literature on SOC phenomena. 

If the physical quantity Zij^k at each cellular node is defined in terms of an energy 
e, the discretized change of the quantity dz / dt corresponds then to a quantized amount 
of dissipated energy < Ae > per node. The instantaneous energy dissipation rate 
de/dt =< Ae > dV/dt of the system scales then with the instantaneous volume 
change dV{t)/dt, while the total amount of dissipated energy during an avalanche 
event can be computed from the time integral, E =< Ae > J* V{t)dt, assuming that 
a mean energy quantum < Ae > is dissipated per unstable node. 

In astrophysical applications, the observed quantity is usually a photon flux, and 
thus the energy quantity e = hv corresponds to the photon energy at frequency v, 
while < Ae >= Nhi^/At is the photon flux of N photons that are radiated in a vol- 
ume cell Ay during a time step At for a given radiation process. The conversion of the 
intrinsic energy dissipation rate de/dt =< Ae > (dV/dt) (e.g., of magnetic energy 
Eb in a magnetic reconnection process) has then to be related to the amount of emitted 
photons < Ae >= Nhv/ At by a physical model of the dissipation process. The phys- 
ical model of the radiation process may involve a nonlinear scaling between the amount 
of dissipated energy and the number of emitted photons in a specific wavelength range, 
which changes the powerlaw slope of the observed distributions of photon fluxes, com- 
pared with the predictions of the generic quantity -2i,j,fe used in SOC cellular automaton 
simulations. 
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2.2 Analytical SOC Models 

The complexity of the spatio-temporal pattern of a SOC avalanche is produced by a 
simple mathematical re-distribution rule that defines the next-neighbor interactions on 
a microscopic level, while the observed structure manifests itself at the macroscopic 
level. While cellular automaton models operate on the microscopic level of a lattice 
cell, the complex macroscopic structure cannot be analytically derived from the micro- 
scopic states of the nonlinear system. In classical thermodynamics, the macroscopic 
state of a gas, such as the velocity distribution function of molecules can be derived 
from the (binomial) probability distribution function of microscopic states. In con- 
trast, nonlinear dissipative systems in the state of self-organized criticality seem to 
exhibit a higher level of complexity, so that their macroscopic morphology cannot be 
analytically derived from the microscopic states. Therefore, analytical models can de- 
scribe the macroscopic SOC parameters only by simplified approximations of the SOC 
dynamics, which are aimed to be consistent with the observed (powerlaw-Uke) proba- 
bility or occurrence frequency distributions. We will describe a few of such analytical 
models in the following, which make quantitative predictions of the powerlaw slopes 
and scaling laws between SOC parameters. Such analytical models can then be used 
for Monte-Carlo simulations of SOC avalanches and be forward-fitted to the observed 
distributions of SOC parameters. 

2.2.1 Exponential-Growth SOC Model (EG-SOC) 

SOC avalanches have always a growth phase with nonhnear characteristics, similar 
to a multiplicative chain-reaction. Mathematically, the simplest function that grows 
in a multiplicative way is the exponential function with a positive growth rate. It is 
therefore not surprising that the first analytical models of SOC processes were based 
on such an exponential growth function (even when they were not called SOC models 
at the pre-Bak time). The attractive feature for SOC applications is the fact that the 
simple assumption of an exponential growth curve combined with random durations 
automatically leads to a powerlaw distribution function, as we will see in the following. 

The earliest analytical models in terms of an exponential growth phase with satu- 
ration after a random time interval go back to WilUs and Yule (1922) who applied it 
to geographical distributions of plants and animals. Yule's model was applied to cos- 
mic rays (Fermi 1949), to cosmic transients and solar flares (Rosner and Vaiana 1978; 
Aschwanden et al. 1998), to the growth dynamics of the world-wide web (Huberman 
and Adamic 1999), as well as to the distribution of the sizes of incomes, cities, inter- 
net files, biological taxa, and in gene family and protein family frequencies (Reed and 
Hughes 2002). 

Here we describe the analytical derivation of the exponential-growth SOC model, 
following Aschwanden (2011a, Section 3.1). We define the time evolution of the en- 
ergy release rate W{t) of a nonlinear process that starts at a threshold energy of Wq 
by 

W{t) = Wo exp ^-^^ , 0<t<T, (2.6) 
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Figure 2.3: Time evolution of energy release rate W{t) for 3 different ratios of growth times 
to saturation times, ra/ts = (0.5,1.0,2.0) (left) and the corresponding powerlaw distribu- 
tions of the peak energy release rate P. Note that the event set with the shortest growth time 
(ra/ts = 0.5) reaches the highest energies and thus produces the flattest powerlaw slope 
{a = l + TG/ts = 1.5}. 



where tq represents the exponential growth time. The process grows exponentially 
until it saturates at time t = t with a saturation energy Ws, 



Ws = W{t = t) = Wo exp — 



(2.7) 



We define a peak energy release rate P that represents the maximum energy release rate 
Ws, after subtraction of the threshold energy Wq, that corresponds to the steady-state 
energy level before the nonhnear growth phase, 



P = Ws-Wo = Wo 





- 1 









(2.8 



In the following, we will refer to the peak energy release rate P also briefly as "peak 
energy". For the saturation times r, which we also call "rise times", we assume a 
random probabiUty distribution, approximated by an exponential function N{t) with 
e-folding time constant ts, 



j>i[T)aT = — exp 
ts 



ts) 



(2.9) 



This probability distribution is normalized to the total number of A^o events. 

In order to derive the probabihty distribution N{P) of peak energy release rates P, 
we have to substitute the variable of the peak energy, P, into the function of the rise 
time t(P), which yields (using the functional relationship t(P) from Eq. 2.8), 



N{P)dP = N{T)dT = N[t{P)] 



dP 



dP = 



No{ap - 1) / P 



Wo 



Wo 



+ 1 



dP , 
(2.10) 
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Figure 2.4: Schematic of the time evolution of an avalanche event, consisting of (i) a rise 
time (r) with exponential growth of the energy release W(t) from a threshold level Wo to the 
saturation level Ws, and (ii) a decay time (D) with a constant decay rate = dW/di = Wo/td. 

which is an exact powerlaw distribution for large peak energies (P Wq) with a 
powerlaw slope ap of 



The powerlaw slope thus depends on the ratio of the growth time to the e-folding sat- 
uration time, which is essentially the average number of growth times. Examples of 
time series with avalanches of different growth times (tg /tg = 0.5, 1 .0, 2.0) are shown 
in Fig. 2.3, along with the corresponding powerlaw distributions of peak energies P. 
Note that the fastest growing events produce the flattest powerlaw distribution of peak 
energies. 

Once an instability has released a maximum amount Ws of energy, say when an 
avalanche reaches its largest velocity on a sandpile, the energy release gradually slows 
down until the avalanche comes to rest. For sake of simplicity we assume a linear decay 
phase of the released energy (Fig. 2.4), 



where t2 is the end time of the process at t2 = ti + D. The time interval T of the total 
duration of the avalanche process is then the sum of the exponential rise phase r and 
the linear decay phase D as illustrated in Fig. 2.4, 



We see that this relationship predicts an approximate proportionaUty of T oc P for 
large avalanches, since the second term, which is linear to P, becomes far greater than 
the first term with a logarithmic dependence (oc In P). 




(2.11) 




(2.12) 




(2.13) 



2.2. ANALYTICAL SOC MODELS 



11 



t= 1 



t=2 



t=3 



■■■■■■■■■ 


■■ 


■■■■■■■■■ 


■■ 


■■■■■ ■■■ 


■■ 


■■■■ ■■ 


■■ 


■■■fe MB 


■■ 




■■■■ 




■■■■■ ■■■ 


■■ 


■■■■■■■■I 


■■ 


■■■■■■■■■ 


■■ 


■■■■■■■■I 


■■ 



Ik jm 



t=4 



t=5 



1^=1 



a. 



U=2 



ts=3 




ts=4 




ts=5 



Figure 2.5: Spatial patterns of a propagating avalanche in subsequent time steps in a 2-D cel- 
lular automaton model with a next-neighbor redistribution rule (top) and time profiles of energy 
release rate (bottom), for saturation times of ts = 1, 2, 5At. The black cells represent cells 
with random fluctuations below the threshold, Zk < Zc, the grey cells contain possibly unstable 
cells with fluctuation z^ > Zc that are subject to a first redistribution, while the white cells have 
already been affected by a redistribution rule before. Most avalanches die out after step t > 2. 



For the calculation of the distribution N{t) we express the total duration T in terms 
of the rise time r and find a powerlaw function for the distribution of flare durations T, 



N{T)dT = N[t{T)] 



dT 



dT 



dT. 



Nojar - 1) /T_ 



1 



TD 



(2.14) 



We define also the total released energy E by the time integral of the energy release 
rate W{t) during the event duration T, but neglect the rise time r (i.e., T k. D) and 
subtract the threshold level Wo before the avalanche, 



E 



f[W{t) 

JO 



Wo] dt 



T+D 



[Wit) 



Wo] dt = -PD . 



(2.15) 



leading to a powerlaw-Uke function for the frequency distribution of energies E, 



N{E)dE = N[P{E)] 



dP 



dE 



dE 



Nojap - 1) 
2Ea 







— OLp 


■ E ' 


V Eo 






_Eo_ 



-1/2 



(2.16) 

Thus, we find the foUowing approximate scaling laws between the powerlaw indices, 

ap = l + Tc/ts 

ar = ap (2-17) 

as = (ap + l)/2 

In summary, this model predicts powerlaw distribution functions for the three SOC 
parameters P, E, and T, which match the simulated distributions with cellular automa- 
ton simulations, as well as the observed distributions of solar flare hard X-ray fluxes 
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(Lu and Hamilton 1991). For a ratio of Ta/ts = 1, this model predicts ap = 2.0, 
ar = 2.0, and as = 1-5. This particular time ratio of Ta/ts = 1 imphes that an 
avalanche typically saturates after one exponential growth time (see the spatio-temporal 
patterns for small avalanches in Fig. 2.5). Open questions are: Which physical pro- 
cess would explain this particular time scale ratio rc/ts = 1 ? Why does the decay 
phase has a Unear behavior? Is the exponential distribution of avalanche growth times 
consistent with observations, since powerlaw-like distributions are expected for SOC 
parameters? How can the observed intermittent fluctuations of time profiles and the 
fractal geometry be accomodated in a model with a monotonic growth function? Al- 
though this model seems not to reproduce all observed properties of SOC phenomena, 
it has a didactical value, since it represents the most basic model that links the nonlinear 
evolution of instabilities to the powerlaw distributions observed in SOC phenomena. 

The exponential-growth model is most suitable for multiplicative avalanche pro- 
cesses, where the increase per time step during the rise phase is based on a multi- 
plicative factor, such as it occurs in nuclear chain reactions, population growth, or ur- 
ban growth. Alternatively, avalanche processes that continuously expand in space may 
show an energy increase that scales with the area or volume, which has a powerlaw rela- 
tionship to the spatial or temporal scale, i.e., A{t) oc r{t)^ oc or V{t) oc r{t)^ oc t^. 
Such a model with a powerlaw-growth function 



rather than an exponential-growth function (Eq. 2.6) was computed in Aschwanden 

(2011a; Section 3.2), which predicts identical scaling laws between the SOC parame- 
ters (P, E, T), but the occurrence frequency distributions exhibit an exponential fall-off 
at the upper end. Otherwise it has similar caveats as the exponential- growth model (i.e., 
monotonic growth curve rather than intermittency, and Euclidean rather than fractal 
avalanche volume). 

Another variant of the exponential-growth model is the logistic-growth model (As- 
chwanden 2011a; Section 3.3), which has a smoother transition from the initially ex- 
ponential growth to the saturation phase, following the so-called logistic equation, 



where the dissipated energy is limited by the so-called carrying capacity limit E^c used 
in ecological applications. The resulting occurrence frequency distributions are pow- 
erlaws for the energy E and peak energy dissipation rate P, but exponential functions 
for the rise time r and duration T of an avalanche. Like the exponential-growth and 
the powerlaw-growth model, the temporal intermittency and the geometric fractality 
observed in real SOC phenomena are not reproduced by the smoothly-varying time 
evolution of the logistic-growth model. 




(2.18) 



dE{t) _ E{t) \ E{t) 



(2.19) 



dt TG _ Eoo 
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2.2.2 The Fractal-Diffusive SOC Model (FD-SOC) 



The microscopic structure of a SOC avalanche has been simulated with a discretized 
mathematical re-distribution rule, which leads to highly inhomogeneous, filamentary, 
and fragmented topologies during the evolution of an avalanche. Therefore, it appears 
to be adequate to develop an analytical SOC model that approximates the inhomoge- 
neous topology of an avalanche with a fractal geometry. Bak and Chen (1989) wrote 
a paper entitled "The physics of fractals", which is summarized in their abstract with 
a single sentence: Fractals in nature originate from self-organized critical dynamical 
processes. 

A statistical fractal-diffusive avalanche model of a slowly-driven SOC system has 
been derived in Aschwanden (2012). This analytical model represents a universal 
(physics-free) description of the statistical time evolution and occurrence frequency 
distribution function of SOC processes. It is based on four fundamental assumptions: 
(1) A SOC avalanche grows spatially like a diffusive process; (2) The spatial volume 
of the instantaneous energy dissipation rate is fractal; (3) The time-averaged fractal 
dimension is the mean of the minimum dimension Ds,min ~ 1 (for a sparse SOC 
avalanche) and the maximum dimension Ds,max = S (given by the Euclidean space); 
and (4) The occurrence frequency distribution of length scales is reciprocal to the size 
L of spatial scales, i.e., N{L) oc in Euclidean space with dimension S. We will 
discuss these assumptions in more detail in the following. 

The first assumption of a diffusive process is based on numerical simulations of cel- 
lular automaton models. A SOC avalanche propagates in a cellular automaton model 
by next-neighbor interactions in a critical state, where energy dissipation propagates 
only to the next neighbor cells (in a S-dimensional lattice grid) that are above a critical 
threshold. This mathematical rule that describes the entire dynamics and evolution of 
a SOC avalanche is very simple for a single time step, but leads to extremely com- 
plex spatial pattems after a finite number of time steps. For a visualization of a large 
number of such complex spatial pattems generated by a simple iterative mathematical 
redistribution rule see, for instance, the book "A New Kind of Science" by Wolfram 
(2002). The complexity of these spatial patterns can fortunately be characterized with 
a single number, the fractal dimension Ds. If one monitors the time evolution of a spa- 
tial pattem of a SOC avalanche in a cellular automaton model, one finds that the length 
scale x{t) evolves with time approximately with a diffusive scaling (see radius r{t) 
of snapshots of a 2-D cellular automaton evolution in Fig. (2.6) and its time evolution 
r{t) oc i^/^ in Fig. (2.8), bottom right panel). 



which leads to a statistical scaling law between the avalanche sizes L = x{t = T) and 
time durations T of SOC avalanches. 



The second assumption of a fractal pattern of the instantaneous energy rate is also 
based on tests with cellular automaton simulations (see measured fractal dimensions 
D2 of snapshots in Fig. (2.7) and the time evolution D2 {t) in Fig. (2.8) top right panel). 



x{t) oc t^'^ 



(2.20) 



L oc . 



(2.21) 
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Time t=26 



Time t=118 



Time t=132 




Figure 2.6: Time evolution of a large avalanche event in a 2-D cellular automaton 
simulation with grid size N = 64^. The 12 panels show snapshots at particular burst 
times from t = 26 to t = 690 when the energy dissipation rate peaked. Active nodes 
where energy dissipation occurs at time t are visualized with black and grey points, 
depending on the energy dissipation level. The starting point of the avalanche occurred 
at pixel (x, y) = (41, 4), which is marked with a cross. The time-integrated envelop 
of the avalanche is indicated with a solid contour, and the diffusive avalanche radius 
r{t) = t^l"^ is indicated with a dashed circle (Aschwanden 2012). 



The fractal dimension is essentially a simplified parameter that describes the "micro- 
roughness", "graininess", or inhomogeneity of critical nodes in a lattice grid in the 
state of self-organized criticality. Of course, such a single number is a gross over- 
simplification of a complex system with a large number of degrees of freedom, but 
the numerical simulations confirm that avalanche pattems are fractal (Fig. 2.7). Thus 
we define the volume Vs(t) of the instantaneous energy dissipation rate in terms of a 
fractal (Hausdorff) dimension Dg that scales with the length scale x as, 

Fs(t) oc ^ (2.22) 



which leads also to a statistical scaling law between avalanche volumes V and spatial 
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Figure 2.7: Determination of the fractal dimension D2 = \ogAi/logXi for the in- 
stantaneous avalanche sizes of 6 time steps of the avalanche event shown in Fig. (2.6). 
Each row is a different time step and each column represents a different binning of 
macropixels (Axi = 1,2,4, 8). The fractal dimension is determined by a linear regres- 
sion fit shown on the right-hand side. The mean fractal dimension of the 12 avalanche 
snapshots shown in Fig. (2.6) is D2 = 1.43 ± 0.17 (Aschwanden 2012). 
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Figure 2.8: Time evolution of the same large avalanche event from a 2-D cellular au- 
tomaton simulation with grid size N = 64^ as shown in Fig. (2.6). The time profiles 
include the instantaneous energy dissipation rate f{t) = de/dt (top left), the time- 
integrated total energy e{t) (bottom left), the instantaneous fractal dimension D2{i) 
(top right), and the radius of the avalanche area r{t) (bottom right). The observed 
time profiles from the simulations are outlined in solid Unestyle and the theoretically 
predicted average evolution in dashed linestyle. The statistically predicted values of 
the instantaneous energy dissipation rate f{t) oc t^/^ (dotted curve) and peak energy 
dissipation rate p{t) cx (dashed curve) after a time interval t are also shown (top left 
panel). The 12 time labels from 26 to 690 (top left frame) correspond to the snapshot 
times shown in Fig. (2.6) (Aschwanden 2012). 



scales L or durations T of SOC avalanches (with Eq. 2.21) 

Vs a oc . (2.23) 

The third assumption of the mean fractal dimension has also been confirmed by 
numerical simulations of cellular automaton SOC processes in all three dimensions 
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S* = 1, 2, 3 (Aschwanden 2012), but it can also be understood by the following plausi- 
bility argument. The sparsest SOC avalanche that propagates by next-neighbor interac- 
tions is the one that spreads only in one spatial dimension, and thus yields an estimate 
of the minimum fractal dimension of Ds.min ~ 1, while the largest SOC avalanche 
is almost space-filling and has a volume that scales with the Euclidean dimension, 
Ds,max = S. Combining these two extremal values, we can estimate a time-averaged 
fractal dimension < Ds > from the arithmetic mean, 

.... Ds ,min + Ds .max ± /n nA\ 

< L>s >~ ^ = , (2.24j 

which yields a mean fractal dimension of < -D3 >= (1 + 3)/2 = 2.0 for the 3D case 
(5 = 3). 

The fourth assumption on the size distribution is a probability argument. The sys- 
tem size Lgys of a SOC system represents an upper Umit of spatial scales L for SOC 
avalanches, i.e., L < Lgys- For the 3D-case, the volumes V of individual avalanches 
are also bound by the volume Vsys of the system size, i.e., V = < Vsys = L^sys- 
If the entire system is in a critical state, SOC avalanches can be produced everywhere 
in the system, and the probability N{L) for a fixed avalanche size L with volume V is 
simply reciprocal to the size, (visuaUzed in Fig. 2.9), i.e., 

N{V) oc oc (2.25) 



V 



which is equivalent to 



N{L) cx^o^L-K (2.26) 

Based on these four model assumptions we can now quantify the time evolution of 

SOC parameters. For sake of simplicity we apply this fractal -diffusive SOC model 
to an astrophysical source that emits a photon flux f{t) that is proportional to the 
instantaneous energy dissipation volume Vs{t) of a SOC avalanche, where a mean 
energy quantum < /^E > is emitted per volume cell element AV. Thus, the emitted 
flux is (using the diffusive scaling x{t) oc t^^'^), 

f{t) = oc< AE > Vs{t) =< AE > x{t)^' =< AE > t^'/'^ . (2.27) 

In the 2-D case with Ds = D2 = {I + 2)/2 = 3/2 (Eq. 2.24) we expect a statistical 
scaUng of f{t) oc t^/^ (see Fig. 2.8 top left). In the 3-D case with Ds = D3 = 
(1 + 3)/2 = 2, we expect than the proportionality f{t) oc t^. In Fig. (2.10) (second 
panel) we simulate such a flux time profile f{t) by applying noise fluctuations in the 
time evolution of the fractal dimension D2{t) (Fig. 2.10, top). 

The statistical peak value p{t) of the energy dissipation rate after time t can be 
estimated from the largest possible avalanches, which have an almost space-filling di- 
mension Ds ^ S, 

= ^ oc< A£; > Vs,maAt) =< Ai? > x{tf =< AE > t^'^ . (2.28) 
at 
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Figure 2.9: Schematic diagram of the Euclidean volume scaling of the diffusive 
avalanche boundaries, visualized as circles or spheres in the three Euclidean space 
dimensions S = 1, 2, 3. The Euclidean length scale x of subcubes decreases by a fac- 
tor 2 in each step (xi = 2~', i = 0, 1, 2), while the number of subcubes increases by 
rii = (2*)'^, defining a probability of N{xi) oc for each avalanche size with size 

Xi . 
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and thus would be expected to scale as p{t) cx t'^/^ cx t^-^ for the S = 3 case. Thus the 
time profile of peak values envelopes the maximum fluctuations of the flux time profile 
f{t) (Fig. 2.10, second panel). 

The evolution of the total dissipated energy e{t) after time t is simply the time 
integral, for which we expect 

e(t) = /* ^lldT oc /* r^-/2 oc ^ (2.29) 

Jo dT Jq 

which yields the function e{t) oc f''/^ for the 2-D case (Fig. 2.8, bottom left), and a 
function e{t) oc t'^ for the 3-D case (Fig. 2.10, third panel). These time evolutions apply 
to every SOC model that has an emission f{t) proportional to the fractal avalanche vol- 
ume Vs (t) . For apphcations to observations in a particular wavelength range there may 
be an additional scaUng law between the avalanche volume and emission (or intensity) 
of the underlying radiation mechanism. 

The fourth assumption on the probability distribution of (avalanche) length scales, 

N{L) oc Vg^ oc L-^ . (2.30) 

is a simple probability argument, implying that the number or occurrence frequency 

of avalanches has an equal likelyhood throughout the system, so it assumes a homo- 
geneous distribution of critical states across the entire system. This is an extremely 
important assumption, which automatically predicts a powerlaw distribution for length 
scales. Using the scaling laws that result from Eq. (2.21) and (2.27)-(2.29) for L = 
x{t = T),E = e{t = T),F = f{t = T), and P = p{t = T) for an avalanch duration 
timeT, 

L oc Ti/2 
F oc 

P oc PS/^ ' (2-^^) 
E a 

we can directly calculate the occurrence frequency distributions for all these SOC 
parameters, by substituting the variables from the correlative relationships given in 
Eq. (2.31), yielding 



N{T)dT = N{L[T\) 



dL 



dT 



dP oc T-[(i+^)/2l dT . 



N{F)dF = N{P[P] 

N{P)dP = N{P[P 
dP 



dP 



dF 



dF oc F-[i+(^-i)/^sl dF , 
dP 



dP 



dP oc p-P-VS] dP , 



N{E)dE = N{P[E]) 



dE 



dE oc £;-[i+(s-i)/(i?s+2)] . 



(2.32) 

(2.33) 
(2.34) 
(2.35) 



This derivation yields naturally powerlaw functions for all parameters L, T, F, P, and 
E, which are the haUmarks of SOC systems. In summary, if we denote the occurrence 
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Table 2.2: Theoretically predicted occurrence frequency distribution powerlaw slopes 

a and power indices /? of parameter correlations predicted for SOC cellular automatons 
with Euclidean space dimensions S = 1,2,3, for the fractal dimension Ds, length 
scale L, time duration T, instantaneous energy dissipation rate (or flux) F, peak energy 

dissipation rale (or peak flux) P. and total time-integrated energy E. 



Theory 


S=l 


S=2 


S=3 


Ds = {l + S)/2 


1 


3/2 


2 


aL = S 


1 


2 


3 


ar = (1 + S)/2 


1 


3/2 


2 


ap = l + {S- 1)1 Ds 


1 


5/3 


2 


ap = 2-1/5' 


1 


3/2 


5/3 


aE = l + {S-l)/{Ds + 2) 


1 


9/7 


3/2 


L oc TV2 


L oc rV2 


L oc TV2 


L oc TV2 


F oc 


F oc 


F oc 


FocTi 


P (X 


P oc Ti/2 


P (xT^ 


p oc r3/2 


E oc 


E oc T3/2 


E oc 


EocT^ 



frequency distribution N{x) of a parameter a; with a powerlaw distribution with power 
index a^, 

N{x)dx oc X""- da; , (2.36) 

we have the following powerlaw coefficients for the parameters x = T,F, P, and 
E, 

ap ={l + S)/2 

ap =1 + {S-1)/Ds , . 

ap =2-1/ S ■ ^^■'^'> 

ap =l + {S-l)/{Ds + 2) 

The powerlaw coefficients and correlation are summarized in Table 2.2 separately 
for each Euclidean dimension S — 1,2,3. 

These correlation coefficients and powerlaw indices of frequency distributions have 
been found to agree within 10% with numerical simulations of cellular automatons 
for all three Euclidean dimensions (5* = 1,2,3) (Aschwanden 2011a). Some devia- 
tions, especially fall-offs at the upper end of powerlaw distributions, are Ukely to be 
caused by finite-size effects of the lattice grid. 



2.2.3 Astrophysical Scaling Laws 

SOC theory applied to astrophysical observations covers many different wavelength 
regimes, for instance gamma-rays, hard X-rays, soft X-rays, and extreme ultra-violet 
(EUV) in the case of solar flares. A comprehensive review of such studies is given in 
Section 7 of Aschwanden (201 la). However, since each wavelength range represents 
a different physical radiation mechanism, we have to combine now the physics of the 
observables with the (physics-free) SOC statistics. 
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Figure 2. 10: Simulation of the fractal-diffusive SOC model for an Euclidean dimension 
S = 3, showing the time evolution of the fractal dimension Ds{t) (top panel), the in- 
stantaneous energy dissipation rate f{t) and peak energy dissipation rate p{t) (second 
panel), the total time-integrated dissipated energy e{t) (third panel), and the soft X-ray 
time profile fsxr{t) (bottom panel), which results from the convolution of the instanta- 
neous energy dissipation rate f{t) (second panel) with an exponential decay function 
with an e-folding time of Tdecay (shown in insert of bottom panel) (Aschwanden and 
Freeland 2012). 



Let us consider soft and hard X-ray emission in solar or stellar flares. Soft X- 
ray emission during solar flares is generally believed to result from thermal free-free 
and free-bound radiation of plasma that is heated in the chromosphere by precipita- 
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fsxr{t) = / fit) CXp 
J —oo 



dt' , (2.38) 



tion of non-thermal electrons and ions, and which subsequently flows up into coronal 
flare (or post-flare) loops, a process called "chromospheric evaporation process" (for 
a review see, e.g., Aschwanden 2004). Therefore, we can consider the flare-driven 
chromospheric heating rate as the instantaneous energy dissipation process of a SOC 
avalanche, as shown in the simulated function f{t) in Fig. 2.10 (second panel). The 
heated plasma, while it fills the coronal flare loops, loses energy by thermal conduction 
and by radiation of soft X-ray and EUV photons, which generally can be characterized 
by an exponential decay function after an impulsive heating spike. In Fig. 2.10 (bot- 
tom) we mimic such a soft X-ray radiation Ught curve by convolving the instantaneous 
energy dissipation rate f{t) (Fig. 2.10, second panel) with an exponentially decaying 
radiation function (with an e-folding time constant of Tdecay), 

'^decay 

which shows also a time dependence that follows approximately 

/.xr(t)oc/(t)octi-V (2.39) 

because the convolution with an exponential function with a constant e-folding time 
constant acts like a constant multiplier In the limit of infinitely long decay times 
(Tdecay oo), OUT convolution function (Eq. 2.38) turns into a time integral of the 
heating function f{t), which is also known as Neupert effect (Dennis and Zarro 1993; 
Dennis et al. 2003). 

The heating function is identified with the non-thermal hard X-ray emission, which 
indeed exhibits a highly fluctuating and intermittent time profile for energies of > 25 
keV, where non-thermal emission dominates, 

fhxr{t)(xp{t)^t'-' , (2.40) 

Thus the occurrence frequency distributions of fluxes are expected to be different for 
soft X-rays and hard X-rays, the hard X-ray flux fhxr{t) follows the statistics of the 
highly fluctuating peak energy dissipation rate p{t), while the soft X-ray flux fsxrit) 
is expected to follow the statistics of the smoothly-varying (convolved) time profile 
f{t). The total duration T of energy release of an avalanche (or flare here) corresponds 
essentially to the rise time trise of the soft X-ray flux, or to the total flare duration for 
hard X-rays, because the decay phase of a soft X-ray flare light curve is dominated by 
conductive and radiative loss, rather than by continued heating input. Thus based on 
the generic relationships sunmiarized in Table 2.2 we expect for the 3-D case (S = 3), 

N{T) oc T-"-^ = , (2.41) 

N{fsxr) oc F-°"' = i?'-[i+(S-i)/i?s] = p-2 ^ (2.42) 

Nifhxr) oc P-""- = p-P-i/s] = p-5/3 ^ (2.43) 

Applications to observations did show satisfactory agreement with these theoretical 
values, i.e. powerlaw slopes of ap = 2.0 for soft X-rays and ap = 1.67 for hard X- 
rays (Aschwanden 201 la; 201 lb; Aschwanden and Freeland 2012), except for the oc- 
currence frequency distributions of flare durations T during solar cycle maxima, when 
the flare pile-up bias appears to have a steepening side-effect (Aschwanden 2012). 
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2.2.4 Earthquake Scaling Laws 

While the foregoing discussion is relevant to astrophysical SOC phenomena (solar and 
stellar flares), similar physical scaling relationships between the observables and SOC 
cellular automaton quantities Zij^k can be developed in other fields. For earthquakes 
in geophysics for instance, measured quantities include the length Ls and width 
of a rupture area, so the rupture area has the scaling A oc LsLyj. For large ruptures, 
however, when the surface rupture length is much larger than the rupture width, i.e., 
Ls ^ iu>, the width L^, was found to be approximately constant, so that large rup- 
tures saturate the faults surface and scale with a fractal dimension of -D = 1, whereas 
a smaller rupture propagates across both dimension {D = 2) of the fault face (Yo- 
der et al. 2012). The modeling of size distributions of earthquake magnitudes can 
then be conducted with a similar methodology as outUned with our fractal-diffusive 
model (Section 2.2.2), requiring: (i) The definition of the Euclidean space dimension 
for earthquakes, which could be 5 = 2 if they are treated as surface phenomena without 
significant depth variabihty, or by 5 = 3 otherwise; (ii) The measurement of the fractal 
dimension D, which can exhibit multi-fractal scaling from £> = 1 for large earthquakes 
to _D = 2 for small earthquakes; and (iii) assigning the observable m, which is the mag- 
nitude of the earthquake, to the corresponding physical quantity, i.e., the peak energy 
dissipation rate P, the smoothed energy dissipation rate F, or total time-integrated en- 
ergy E. Furthermore, earthquake size distributions are often reported in terms of a 
cumulative distribution function A^(> m), which has a powerlaw index /3 « a — 1 that 
is flatter by one than the differential distribution function N{m)dm. 

2.3 Alternative Models Related to SOC 

Here we discuss a number of alternative dynamical models that are related to SOC 
models or have similar scaling laws, and discuss what they have in common with SOC 
or where they differ. A metrics between processes and observables is synthesized in 
Fig. 2.15. 

2.3.1 Self-Organization Without Criticality (SO) 

Self-organization (SO) is often referred to geometric patterns that originate by mu- 
tual interaction of their elements, without coordination from outside the system. For 
instance, convection arranges itself into a regular pattern of almost equal-sized con- 
vection cells, a structuring process also known as Benard cells. Other examples are 
the regular wavy pattern of sand dunes in the desert, wavy pattems of Cirrus clouds in 
the Earth atmosphere, Jupiter's atmosphere with white bands of ammonia ice clouds, 
the spiral-like pattems of the Belousov-Zhabotinsky reaction-diffusion system, or ge- 
ometric patterns in biology, such as the skin of zebras, giraffes, tigers, tropical fishes, 
or formation flight of birds. From these examples we see that self-organization mostly 
refers to the self-assembly of geometric pattems, which are more or less stable over 
long time intervals, although they arise from the physics of non-equilibrium processes, 
which can involve diffusion, turbulence, convection, or magneto-convection, which are 



24 



CHAPTER 2. THEORETICAL MODELS OF SOC SYSTEMS 



governed by long-range interactions (via pressure and forces). 

What is the difference to self-organized criticality (SOC)? Are sandpile avalanches 
a self-organizing (SO) pattern? In the standard scenario of the sandpile SOC model, 
individual avalanches are a local phenomenon that are randomly triggered in space and 
time, but occur independently, at least in the slowly-driven case. Thus, one avalanche 
has no mutual interaction with another avalanche and the outcome of the final size 
is independent of another, in contrast to self-organization without criticality, where 
the interaction between system- wide structures is coupled. In other words, sandpile 
avalanches are govemed by localized disturbances via next-neighbor interactions, while 
self-organizing pattems may be formed by both next-neighbor and long-range interac- 
tions. Another difference is that SO creates spatial patterns, while SOC generates dy- 
namical events (i.e., avalanches). Also the statistical distributions of the two processes 
are different. A self-organizing pattern is likely to produce a preferred size scale (e.g., 
the solar granulation or the width of zebra stripes), while self-organized criticality pro- 
duces a scale-free powerlaw distribution of avalanche sizes. The difference between 
SO and SOC could may be best illustrated with a sandpile analogy. The same sandpile 
can be subject to self-organization (SO), for instance when a steady wind blows over 
the surface and produces wavy ripples with a regular spacing pattern, as well as be 
subject to self-organized criticality (SOC), when intermittent avalanches occur due to 
random-like disturbances by infalling sand grains. The former geometric pattem may 
appear as a spatial-periodic pattem, while the latter may exhibit a fractal geometry. 

2.3.2 Forced Self-Organized Criticality (FSOC) 

In Bak's original SOC model, which is most genuinely reproduced by sandpile avalan- 
ches and cellular automaton simulations, criticality is continuously restored by next- 
neighbor interactions. If the local slope becomes too step, an avalanche will erode it 
back to the critical value. If the erosion by an avalanche flattened the slope too much, it 
will be gradually restored by random input of dropped sand grains in the slowly-driven 
Umit. When the observations were extended to magnetospheric substorms in the night- 
side geotail, powerlaw-like size distributions were found, suggesting a self-organized 
criticality process, but at the same time correlations with magnetic reconnection events 
at the dayside of the magnetosphere were identified, driven by solar wind fluctuations, 
suggesting some large-scale transport processes and long-range coupling between the 
day-side and night-side of the magnetosphere. Thus, these long-range interactions that 
trigger local instabilities in the current sheet in the Earth's magnetotail were interpreted 
as external forcing or loading, and a combined forced and/or self-organized criticality 
(FSOC) process was suggested (Chang 1992, 1998a,b, 1999a,b). A forced SOC pro- 
cess may apply to a variety of SOC phenomena in the magnetosphere, such as mag- 
netotail current disruptions (Lui et al. 2000), substorm current disruptions (Consolini 
and Lui 1999), bursty bulk flow events (Angelopoulos et al. (1996, 1999), magnetotail 
magnetic field fluctuations (Hoshino et al. 1994), auroral UV blobs (Lui et al. 2000; 
Uritsky et al. 2002), auroral optical blobs (Kozelov et al. 2004), auroral electron (AE) 
jets (Takalo et al. 1993; Consolini 1997), or outer radiation belt electron events (Crosby 
et al. 2005). 

What distinguishes the FSOC model from the standard (BTW-type) SOC model 
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Figure 2.11: Numerical simulation of ID Brownian random walk (left side) and size 
increase of a 2D cellular automaton avalanche (right side) as a function of time t. Both 
processes foUow the trend x{t) oc f^/^ expected for classical diffusion. 




is mostly the long-range coupling, which is absent in sandpiles and cellular automa- 
ton models. The driving force in the FSOC model is a long-distance action (loading 
process), while it is a random local disturbance in the BTW model. Otherwise, the 
FSOC model entails also powerlaw distributions for the avalanche events, fractality, 
intermittency, statistical independence of events (indicated by random waiting time 
distributions), and a critical threshold (for the onset of local plasma instabiUties and/or 
magnetic recormection processes). 

2.3.3 Brownian Motion and Classical Diffusion 

Brownian motion is a concept from classical physics that describes the random mo- 
tion of atoms or molecules in a gas (named after the Scottish botanist Robert Brown), 
which can be best observed when a neutral gas of a different color is released in the 
atmosphere. What we observe then is an isotropic diffusion process (in the absence of 
forces or flows) that monotonically increases with the square-root of time, 

< x{t) >oc t^/^ . (2.44) 

This statistical trend was derived in classical thermodynamics, assuming a Gaussian 
distribution of velocities for the gas molecules, so that the kinetic energy of particles 
follow a Boltzmann distribution in thermodynamic equilibrium. The classical diffu- 
sion process can also be described by a differential equation for a distribution function 
f{x, t) of particles, 

which can describe heat transport, diffusion of gases, or magnetic diffusivity on solar 
and steUar surfaces (i.e., manifested as meridional flows during a solar/steUar activity 
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An example of such a diffusive random walk is simulated in Fig. 2.11 (left panel) 
for the 1-D variable x{t). We have to be aware that a diffusive random walk x{t) of 
a single particle can show large deviations from the expected trend oc t^^^, which is 
only an expectation value for the statistical mean of many random walks < x{t) >= 

Xi{t). In Fig. 11 (right panel) we show also the time evolution of the mean radius 
< r{t) > of a simulated cellular automaton avalanche (taken from Fig. 2.8, bottom 
right), which shows the same trend of a time-dependence of r{t) oc t^^'^. Apparently, 
the enveloping volume of unstable cells in a SOC avalanche, defined by the a math- 
ematical re-distribution rule applied to a coarse-grained lattice with a rough surface 
of stable, meta-stable, and unstable nodes produces a similar time evolution as a ran- 
dom walk (as visuaUzed in Fig. 2.6), leading us to the fractal-diffusive SOC model 
described in Section 2.2.2. This behavior of a diffusive transport process is therefore 
common to both the Brownian motion (or classical diffusion) and to a cellular automa- 
ton SOC process, but it operates for a SOC process only during a finite time interval, 
as long as avalanche propagation is enabled by finding unstable next neighbors, while 
a classical diffusion process goes on forever without stopping. Therefore, we cannot 
define an event and occurrence frequency distributions. However, we can measure the 
fractal dimension of a random walk spatial pattern, which has a Hausdorff dimension 
of 1)2 = 2 in 2-D Euclidean space, or the power spectrum, i.e., P{i^) oc also 
called Brownian noise. 

Generalizations that include power spectra P(i^) oc with arbitrary powerlaw 
indices /3 have been dubbed fractional Brownian motion (fBM), for instance white 
noise (^ = 0, where subsequent steps are uncorrelated), 1/f-noise, flicker noise, or 
pink noise (/? = 1), Brownian noise or red noise ((3 = 2), or black noise (/3 = 3, where 
subsequent time steps have some strong correlations, producing slowly-varying time 
profiles). The latter fluctuation spectrum has been found to describe the stock market 
well. Examples for these different types of fractional Brownian motion are given in 
Fig. 2.12, while the nomenclature of noise spectra is summarized in Table 2.3. 

2.3.4 Hyper-Diffusion and Levy Flight 

The cellular automaton mechanism, the prototype of SOC dynamics, involves a mathe- 
matical re-distribution rule amongst the next neighbor cells (Eq. 2.1). This discretized 
numerical re-distribution rule has been transformed in the continuum hmit to an analyt- 
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Figure 2.12: Noise power spectra (left panels) and corresponding time series (right panels) 
for power spectral indices p = (top row: white noise spectrum), p = 1 (second row: pink 
noise spectrum), p = 2 (third row: red noise spectrum), and p = 3 (bottom row: black noise 
spectrum). The white noise spectrum is multiplied with in the other cases. The time series 
are reconstructed with the inverse Fast Fourier Transform ( Aschwanden 201 la). 



ical function A{x, t), which can be expressed as fourth-order hyper-diffusion equation 
(Liu et al. 2002; Charbonneau et al. 2001), 

where A represents the placeholder of the physical quantity corresponding to the sym- 
bol Zij^k in Eq. 2.1, x is the S-dimensional space coordinate (xij^k for S = 3), 
with second-order centered differencing in space, forward differencing in time, and 
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K = 1/25 the hyper-diffusion coefficient for EucHdean dimension S. However, in 
contrast to classical physics, the hyper-diffusion coefficient k is subject to a threshold 
value in the SOC model, 

« = ( 'lu ^ • (2-46) 

1^ otherwise ^ ' 

So, on a basic level, the evolution of a cellular automaton avalanche translates into 
a hyper-diffusion process, as demonstrated in Liu et al. (2002). How can it be de- 
scribed by classical diffusion as we discussed in the foregoing Section and in the 
fractal-diffusive SOC model in Section 2.2.2? The major difference of the two appar- 
ently contradicting descriptions lies in the fractality: The fourth-order hyper-diffusion 
system (subjected to a threshold instability in the diffusion coefficient k) is applied 
to a EucUdean space with dimension S, while the fractal-diffusive SOC model has a 
second-order (classical) diffusion coefficient in a fractal volume Vs oc x^^ with a frac- 
tal dimension Ds « (1 + 5')/2. The two diffusion descriptions apparently produce an 
equivalent statistics of avalanche volumes Vs, as demonstrated by the numerical simu- 
lations of powerlaw distribution functions N{Vs) by both models (Liu et al. 2002). 

Other modifications of random walk or diffusion processes have been defined in 
terms of the probability distribution function of step sizes. Classical diffusion has 
a normal (Gaussian) distribution of step sizes, which was coined Rayleigh flight by 
Benoit Mandelbrot, while Levy flight was used for a heavy-tailed probabihty distri- 
bution function (after the French mathematician Paul Pierre Levy). Related are also 
the heavy-tailed Pareto probability distribution functions. Thus, Levy flight processes 
include occasional large-step fluctuations on top of classical diffusion, which is found 
in earthquake data, financial mathematics, cryptography, signal analysis, astrophysics, 
biophysics, and solid state physics. A number of SOC phenomena have been analyzed 
in terms of Levy flight processes, such as rice piles (Boguna and Corral 1997), ran- 
dom walks in fractal environments (Hopcraft et al. 1999; Isliker and Vlahos 2003), 
solar flare waiting time distributions (Lepreti et al. 2001), or extreme fluctuations in 
the solar wind (Moloney and Davidsen 2010, 201 1). 

2.3.5 Nonextensive Tsallis Entropy 

Related to Levy flight is also the nonextensive Tsallis entropy, which originates from 
the standard extensive Boltzmann-Gibbs-Shannon (BGS) statistics, where the entropy 
S is defined as, 

S = -kB^Pilnpi, (2.47) 

with ks the Boltzmann constant and pi the probabilities associated with the micro- 
scopic configurations. The standard BGS statistics is called extensive when the cor- 
relations within the system are essentially local (i.e., via next-neighbor interactions), 
while it is called nonextensive in the case when they are non-local or have long-range 
coupling, similar to the local correlations in classical diffusion and non-local steps 
in the Levy flights. As an example, the dynamic complexity of magnetospheric sub- 
storms and solar flares has been described with nonextensive Tsalhs entropy in Balasis 
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et al. (201 1). Although nonlocal interactions are not part of the classical BTW automa- 
ton model, the extensive Tsallis entropy produces statistical distributions with similar 
fractality and intermittency. The nonextensive Tsallis entropy has a variable parameter 
Q, 

^« = ^^(i-Ep?) ' (2.48) 

where g is a measure of the nonextensivity of the system. A value ofq = l corresponds 
to the standard extensive BGS statistic, while a larger value g > 1 quantifies the non- 
extensivity or importance of long-range coupling. A value of g = 1.84 was found to fit 
data of magnetospheric substorms and solar flare soft X-rays (Balasis et al. 201 1). 



2.3.6 Hirbulence 

Turbulence is probably the most debated contender of SOC processes, because of the 
many common observational signatures, such as the (scale-free) powerlaw distributions 
of spatial and temporal scales, the power spectra of time profiles, random waiting time 
distributions, spatial fractality, and temporal intermittency. 

Let us first define turbulence. Turbulence was first defined in fluid dynamics in 
terms of the Navier-Stokes equation. A fluid is laminar at low Reynolds numbers (say 
at Reynolds numbers of i? < 5000, defined by the dimensionless ratio of inertial to 
viscous forces, i.e., R = vL/v^nsc with v the mean velocity of an object relative to 
the fluid, L the characteristic linear dimension of the fluid, and Vyis the kinematic 
viscosity), while it becomes turbulent at high Reynolds numbers. 

SOC phenomena and turbulence have been studied in astrophysical plasmas (such 
as in the solar corona, in the solar wind, or in stellar coronae), where the magneto- 
hydrodynamic (MHD) behavior of the plasma is expressed with the continuity equa- 
tion, the momentum equation and induction equation. 



Dv 



Pg + (j X B) 



iv(V-v) 



(2.49) 



_ = V X (v X B) + r7„ 



y^B 



(2.50) 



with p being the plasma density, p the pressure, B the magnetic field, j the electric 
current density, Vyi^c the kinematic viscosity or shear viscosity, = c? /At:(j the 
magnetic diffusivity, and a the electric conductivity. In the solar corona and in the solar 
wind, the magnetic Reynolds number is sufficiently high {Rm — vL/rjm ~ 10^ — 
10^^) to develop turbulence. Turbulence in the coronal plasma or in the solar wind 
is created by initial large-scale disturbances (for instance by the convective random 
motion in sub-photospheric layers or by the shock-like expansion of flares and coronal 
mass ejections). The large-scale disturbances pump energy into the coronal magnetic 
field or heliospheric solar wind at large scales, which cascade in the case of turbulent 
flow to smaller scales, where the energy can be more efficiently dissipated by friction, 
which is quantified by the kinematic or (Braginskii) shear viscosity coefficient Vyisc- 
Ultimately, the energy of an MHD turbulent cascade is dissipated in the solar wind at 



30 



CHAPTER 2. THEORETICAL MODELS OF SOC SYSTEMS 




lul - ■ — ■ 



0.01 0.10 1.00 10.00 100.00 

Figure 2.13: A spectrum of the solar wind is shown, based on CLUSTER observations from 
large to small scales, with the proton and electron gyroradius scale indicated. The solar wind 
spectrum is interpreted in terms of a turbulent MHD cascade, with the theoretically predicted 
slopes of f~^^^ and f~^^^ from gyro-kinetic theory. The plot proves that the energy continues 
cascading below the proton scale down to the electron scale, where it is converted to heat (via 
electron Landau damping resonance) causing the steepening of the spectrum to Z""* (Howes 
et al. 2008; Sahraoui et al. 2009; credit: ESA, CLUSTER). 

the spatial scale of (gyrating) thermal protons (« 50 km), and at the scale of thermal 
electrons (« 0.5 km). The resulting power spectrum of the solar wind shows a power 
spectrum of P(i') oc u~^/^ at frequencies below the proton scale, P(i') oc 
between the proton and electron scale, and Piu) oc beyond the electron scale 
(Fig. 2.13). For literature references on MHD turbulence in the solar corona and in the 
solar wind see, e.g., Aschwanden (201 la. Section 10.4 therein). 

What is the relationship between turbulence and SOC processes? The transition 
from laminar flow to turbulent flow at a critical Reynolds number Rcrii represents 
a similar thresholded instability criterion as the critical value Zcrit in SOC systems, 
above which a turbulent avalanche could occur with subsequent cascading from vor- 
tices at large scales and little energy dissipation towards smaller scales with stronger 
energy dissipation. Turbulent energy dissipation probably reduces the mean velocity 
of the fluid particles that are faster than the background flow and laminar flow could be 
restored with lower particle velocities and a lower Reynolds number. To make this pro- 
cess self-organizing, we need also a driver mechanism that brings the fluid speed back 
up to the critical Reynolds number. In the solar wind, for instance, there is system- 
atic acceleration with heliocentric distance, which could drive the system from laminar 
back to turbulent flows, and thus it could be self-organizing. Fluctuations of the solar 
wind speed could therefore be considered as SOC avalanches. This could explain the 
scale-free powerlaw distributions of spatial, temporal, and energy scales measured in 
solar wind fluctuations, its fractality and intermittency. If the critical threshold is ex- 
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ceeded only in localized regions, rather than in the entire system in a fully turbulent 
state, energy dissipation would also occur in locally unstable regions, similar to SOC 
avalanches, and thus the two processes may exhibit the same statistical distributions. 
Consequently, turbulence could qualify as a SOC process if it is driven near the laminar- 
turbulent critical Reynolds number in a self-organizing way. However, fully developed 
turbulence, where the entire system is governed by a high (super-critical) Reynolds 
number, would correspond to a fast-driven SOC system in a catastrophic phase with 
permanent avalanching, without restoring the critical state as in a slowly-driven SOC 
system. 

2.3.7 Percolation 

Percolation is a dynamic process that depends on the connectedness and propagation 
probability of next-neighbor elements, and has a lot in common with diffusion, fractal 
structures, and SOC avalanches. A classical paradigm of a percolation process is coffee 
percolation, which contains a solvent (water), a permeable substance (coffee grounds), 
and soluble constituents (aromatic chemicals). The bimodal behavior of a percolation 
process can be best expressed with the observation whether a liquid that is poured on 
top of some porous material, will be able to make its way from hole to hole and reach 
the bottom. The answer will very much depend on the inhomogeneity characteristics 
of the porous material. In the subcritical state, a percolating cluster will exponentially 
die out, while it will propagate all the way to the bottom in a supercritical state. The 
process can be described by branching theory, where the connections between two 
next neighbors can be open and let the liquid pass with probability p, or they can be 
closed and the probability to pass is (1 — p). Completion of a pass from the top to the 
bottom can be expressed as the combined statistical probability along the entire pass. 
The critical value pcrit that decides between the two outcomes is Pcrit = 1/2 for a 
2-D process. So, this process has an extreme sensitivity to the system probability p, 
depending whether it is below or above the critical value. In SOC systems, the outcome 
of an avalanche is independent of the initial conditions of the triggering disturbance. 
Also, the material in a percolation process has a constant value of p and does not self- 
adjust to a critical value, and thus it does not behave like a SOC system. However, 
what a percolation process has in common with a SOC system is the fractality and 
intermittency of propagating features. Other applications of percolation theory are 
used in physics, material science, complex networks, epidemiology, and geography. 

2.3.8 Phase Transitions 

In classical thermodynamics, a phase transition describes the transformation of one 
phase to another state of matter (solid, liquid, gas, plasma). For instance, transitions 
between solid and liquid states are called "freezing" and "melting", between liquid and 
gas are "vaporization" and "condensation", between gas and plasma are "ionization" 
and "deionization". Phase transitions from one state to another can be induced by 
changing the temperature and/or pressure. 

Since there are critical values of temperature and pressure that demarcate the tran- 
sitions, such as the freezing temperature at 0° C or the boihng point of 100° C, we may 
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compare these critical limits with the critical threshold Zcrit in SOC systems. Could 
this critical value be restored in a self-organizing way? To some extent, the daily 
weather changes can be self-organizing. Let us assume some intermediate latitude on 
our planet where the temperature is around the freezing point at some seasonal period, 
say Colorado or Switzerland around November. The temperature may drop below the 
freezing point during night, triggering snowfall, and may raise slightly above freezing 
point during sunny days, triggering melting of snow. In this case, the day-night cycle 
or the Earth's rotation causes the temperature to oscillate around the freezing point, and 
therefore is in a self-organizing state. Additional variation of cloud cover may intro- 
duce fluctuations on even shorter time scales. However, this bimodal behavior of the 
temperature is not a robust operation mode of self-organized criticality, because sea- 
sonal changes will bring the average daily temperature systematically out of the critical 
range around the freezing point. So, such temporary fluctuations around a critical point 
seem to be consistent with a SOC system only on a temporary basis or in intermittent 
(seasonal) time intervals. In fact, a number of metereological phenomena have been 
found to exhibit powerlaw distributions and where associated with SOC behavior, such 
as Nile river fluctuations (Hurst 1951), rainfall (Andrade et al. 1998), cloud formation 
(Nagel and Raschke 1992), climate fluctuations (Grieger 1992), aerosols in the atmo- 
sphere (Kopnin et al. 2004), or forest fires (Kasischke and French 1995, Malamud et 
al. 1998). Phase transitions may be involved in some of these processes (rainfall, cloud 
formation, forest fires). 

Phase transitions are also involved in some laboratory experiments in material or 

solid-state physics, quantum mechanics, and plasma physics, that have been associ- 
ated with SOC behavior, if we extend the term phase transition also to morphological 
changes between highly-ordered and/or chaotic- structured pattems. Such phase tran- 
sitions with SOC behavior have been found in transitions between the ferromagnetic 
and paramagnetic phases, in specimens of ferromagnetic materials settling into one of 
a large number of metastable states that are not necessarily the energetically lowest 
ones (Che and Suhl 1990), in avalanche-like topological rearrangements of cellular do- 
main pattems in magnetic gamet films (Babcock and Westervelt 1990), in the noise 
in the magnetic output of a ferromagnet when the magnetizing force apphed to it is 
charged (Cote and Meisel 1991), called the Barkhausen effect, in superconducting vor- 
tex avalanches in the Bean state (Field et al. 1995), in tokamak plasma confinement 
near marginal stability driven by turbulence with a subcritical resistive pressure gradi- 
ent (Carreras et al. 1996), or in the electrostatic floating potential fluctuations of a DC 
glow charge plasma (Nurujjaman and Sekar-Iyenbgar 2007). 

In summary, many SOC phenomena are observed during phase transitions, pro- 
ducing avalanche-hke rearrangements of geometric patterns, but a phase-transition is 
not necessarily self-organizing in the sense that the critical control parameter stays 
automatically tuned to the critical point. A sandpile maintains its critical slope au- 
tomatically in a slowly-driven operation mode. If a thermostat would be available, 
temperature-controlled phase transitions could be in a self-organized operation mode, 
which indeed can be arranged in many situations. So we conclude that phase transi- 
tions with some kind of "thermostat" tuned to the critical point represent indeed SOC 
systems. 
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Figure 2.14: Network of American airline UNITED with 8 US-based hubs and connections to 
378 destinations in 59 countries (Credit: United). 

2.3.9 Network Systems 

The prototype of a SOC cellular automaton model is a regular lattice grid, where re- 
distributions occur with next neighbors cells. In contrast, networks are irregular nets 
of nodes that are inter-connected in manifold patterns, containing not necessarily only 
next-neighbor connections, but also arbitrary non-local, long-range connections. Pop- 
ular examples of networks are subway maps, city maps, road maps, airplane route 
maps, electric grid maps, social networks, company organizational charts, financial 
networks, etc, all being created by human beings in some way or another. One big 
difference to SOC systems is therefore immediately clear, the addition of long-range 
connections. Furthermore, while the statistical probabilities for next-neighbor interac- 
tions are identical for each lattice point in a SOC cellular automaton model, networks 
can have extremely different connection probabilities at each node point. Nodes that 
have most of the connections are also called hubs. For instance, the American airline 
United has 378 destinations, but server most of them from 10 hubs (Denver, Washing- 
ton DC, San Francisco, Los Angeles, Chicago, Houston, Cleveland, Newark, Narita, 
Guam) (Fig. 14). Since the cormection probability varies greatly from node to node, 
we do expect different clustering patterns in networks and in SOC lattice systems. 

Nevertheless, SOC behavior has been studied in network systems like city growth 
and urban growth (Zipf 1949; Zanette 2007), cotton prizes (Mandelbrot 1963), finan- 
cial stock market (Scheinkman and Woodford 1994; Somette et al. 1996; Feigenbaum 
2003; Bartolozzi et al. 2005), traffic jam (Bak 1996; Nagel and Paczuski 1995), war 
casualties (Richardson 1941, Levy 1983), social networks (Newman et al. 2002), in- 
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ternet traffic (Willinger et al. 2002), and language (Zipf 1949). SOC phenomena with 
powerlaw-like frequency distributions have been observed for phenomena that occur 
in these networks, measured by some quantity that expresses fluctuations above some 
threshold or noise level. Spatial pattems of network phenomena have been found to be 
fractal, and time profiles of these fluctuations were found to be intermittent, all hall- 
marks of SOC systems. So, are those network phenomena self-organizing? What is 
the critical threshold and what drives the system back to the critical threshold? The 
driving force of cotton prices, lottery wins, and stockmarket fluctuations is simply the 
human desire to gain money and to make profit. The driving mechanism for traffic jams 
is certainly the daily need for commuting and transportation in a finite road network. 
The unpredictable number of casualties in crimes and wars is a bit more controver- 
sial, but nobody would disagree that the human desire for power and control plays a 
role. The driving mechanism for social networks is most likely linked to curiosity and 
desire for information. Since all these human traits are quite genuine and persistent, 
a self-organizing behavior is warranted, where the system is continuously driven to- 
wards a critical state, with occasional larger fluctuations (Wall Street crash. Second 
World War, or other "social catastrophes"), on top of the daily fluctuations on smaller 
scales. In summary, most network-related phenomena can exhibit SOC behavior, ap- 
parently with similar powerlaw-like statistics in grids with non-uniform connectivities 
and non-local long-range connections (in network systems) as in regular lattice grids 
with next-neighbor interactions only (in classical SOC systems). 

2.3.10 Chaotic Systems 

Chaotic systems are nonlinear dissipative dynamical processes that are in principle de- 
terministic, but exhibit "chaotic" behavior in the sense of irregular and fractal geometry, 
and intermittent time evolution. Some chaotic systems can be described as simple as 
with two (such as the Lotka-Volterra equation) or three coupled differential equations 
(e.g., the Lorenz equations). Classical examples include forced pendulum, fluids near 
the onset of turbulence, lasers, nonlinear optical devices, Josephson junctions, chemical 
reactions, the three-body system in celestial mechanics, ecological population dynam- 
ics, or the heartbeat of biological organisms (e.g., Schuster 1988). If two dynamical 
variables (say x{t) and y{t)) are plotted in phase space (y versus x), chaotic systems 
often exhibit a cyclic behavior, with a limit cycle that is also called strange attractor. If 
a dynamical system is driven from small disturbances around an equilibrium solution, 
where it exhibits a quasi-linear behavior, bifurcations of system variables x{t) can oc- 
cur at critical values with transitions to chaotic and intermittent behavior (also called 
route to chaos). Such bifurcations ressemble the phase transitions we discussed earlier 
(Section 2.3.8), which we generalized to transitions from highly-ordered to chaotic- 
structured pattems. The complexity of nonlinear chaotic systems is often characterized 
with the dimension of strange attractors (e.g., Grassberger 1985; Guckenheimer 1985), 
which approximately corresponds to the number of coupled differential equations that 
is needed to describe the system dynamics. 

If we compare the dynamics of chaotic systems with SOC processes, a major dif- 
ference is the deterministic evolution of chaotic systems, while SOC avalanches are 
statistically independent of each other, even when they occur subsequently within a 
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small (waiting) time interval. A chaotic system evolves with a system-wide dynamical 
behavior that is influenced by the entire system, as the coupled differential equations 
express it, which includes next neighbors as well as non-local, long-range coupUng, 
while an avalanche in a SOC system occurs only as a result of next-neighbor inter- 
actions. Do chaotic systems have critical thresholds that are self-organizing? Chaotic 
systems do have critical points, where onset of chaos starts, like at pitchfork bifurcation 
point that leads to frequency doubling, but there is usually not a driver that keeps the 
system at this critical point. For instance, deterministic chaos sets in for fluids at the 
transition from the laminar to the turbulent regime, but generally there exists no driv- 
ing mechanism that keeps a chaotic system at the critical Reynolds number, so that the 
system is self-organizing (Section 2.3.6). However, in other respects, chaotic systems 
exhibit similar complexity as SOC systems, regarding fractality and intermittency, and 
even powerlaw distributions may result in the statistics of chaotic fluctuations. For 
weakly nonlinear dissipative systems near the limit cycle, however, small fluctuations 
around a specific spatial scale or temporal scale may produce Gaussian-like distribu- 
tions, e.g., time scales may be centered around the inverse frequency of the limit cycle. 

Let us mention some examples of chaotic behavior in astrophysics. Time series 
analysis of the X-ray variability of the neutron star Her X-1 revealed a low-dimensional 
attractor (Voges et al. 1987; Cannizzo et al. 1990), as well as for the Vela pulsar (Hard- 
ing et al. 1990). Transient chaos was detected for the low-mass X-ray binary star Scor- 
pius X-1 (Scargle et al. 1993; Young and Scargle 1996), as well as for the R scuti star 
(Buchler et al. 1996). The N-body system in celestial mechanics can lead to chaotic 
behavior, such as for the Satumian moon Hyperion (Boyd et al. 1994). Chaotic behav- 
ior with low-dimensional strange attractors and transitions to period doubling was also 
found in solar radio bursts with quasi-periodic time series (Kurths and Herzel 1986, 
Kurths and Karlicky 1989). Chaotic dynamics was found in the solar wind (Polygian- 
nakis and Moussas 1994), in hydrodynamic convection simulations of the solar dynamo 
(Kurths and Brandenburg 1991), and in solar cycle observations (Kremliovsky 1994; 
Charbonneau 2001; Spiegel 2009). Again, all these examples reflect a system-wide 
nonUnear behavior, where subsequent fluctuations are highly correlated in space and 
time, unlike the statistical independency of SOC avalanches. However, chaotic systems 
may exhibit similar statistics of fractal spatial and intermittent time scales at critical 
points and transitions to chaotic dynamics, without having an intrinsic mechanism that 
keeps the chaotic system near this critical point in a self-organizing way. 

2.3.11 Synopsis 

At the end of this chapter we summarize the characteristics of SOC, SOC-related, and 
non-SOC processes in a metrics as shown in Fig. 2.15. The SOC or SOC-like processes 
are listed in the rows of Fig. 2.15, while the observational characteristics are listed in 
the columns. The metrics tabulates whether a specific SOC-like process can exhibit 
the main properties of SOC processes, such as the powerlaw distributions of various 
parameters, the fractal geometry, the temporal intermittency, the statistical indepen- 
dence of events, the restoration of a critical threshold, and next-neighbor interactions. 
The last column contains a property that SOC processes do not have, namely non-local 
and long-range coupUng, which may also represent an important discrimination cri- 
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Fi gure 2.15: Metrics of SOC or SOC-like processes (rows) and observational properties 
(columns), with an evaluation whether a specific process can exhibit a particular observational 
characteristics (Y=Yes, N=No, and the symbol — is filled in for observational properties that do 
not apply). See Chapter 2.3 for details. 



terion between SOC and non-SOC processes, a complementary characteristic to local 
and next-neighbor interaction processes such as the cellular automaton. The properties 
listed in Fig. 2.15 reflect general trends rather than strictly-valid matches. Many pro- 
cesses can exhibit powerlaw distributions of parameters, but there exist always excep- 
tions or deviations from strict powerlaw distributions. The probably most fundamen- 
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tal characteristics of SOC processes is a suitable mechanism that restores the critical 
threshold for a instability, which often is not automatically operating (such as in turbu- 
lence, phase transitions, or chaotic systems), but can be artificially or naturally added 
to a process (such as a "thermostat" for phase transitions). Another important discrim- 
ination criterion between SOC and non-SOC processes is the non-local or long-range 
coupling, which does not exist in the classical SOC cellular automaton mechanism, but 
operates in many other processes (e.g., turbulence, phase transitions, network systems, 
and chaotic systems), while the property of local and next-neighbor interactions exists 
in virtually all processes, and thus is not discriminative. Finally, fractality, intermit- 
tency, and powerlaw behavior is present in most of the processes, so it is not a good 
discrimination criterion between SOC and non-SOC processes, although powerlaws 
have always been considered as the hallmark of SOC processes. 
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